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Chapter  1 


> 


Introduction 


In  this  Report,  we  summarize  results  obtained  with  support  of  the  Air  Force  Contract 
F49620-94-C-0034.  As  a  result  of  this  effort  we  have  developed,  implemented,  and  applied 
several  Eulerian/Lagrangian  computational  models  for  efl&cient  numerical  simulation  of 
high  energy  density  plasmas.  Prototype  setups  studied  in  this  project  are  modeled  after 
physical  systems  encountered  in  advanced  weapons,  space  propulsion  thrusters,  pulsed 
power  systems,  and  other  areas  of  direct  interest  to  the  Air  Force.  The  main  conclusion  of 
our  study  is  that  computationally  efficient  and  physically  sound  description  of  nonsteady 
plasmas  typical  for  these  applications  is  possible  using  the  advanced  hybrid  MHD/PIC 
methods  developed  here. 

In  Section  1,  we  outline  technical  details  of  our  computational  approach,  including 
numerical  algorithms,  physical  models,  and  efficient  numerical  implementation  strategies 
evaluated/implemented  during  this  work.  In  Section  2,  we  discuss  several  representative 
flow  simulations  obtained  using  computational  models  developed  as  a  result  of  our  research. 
In  particular,  we  have  investigated  optimal  regimes  of  plasma  focusing  and  energy  conver¬ 
sion  across  a  range  of  plasma  beam  parameters.  Overall,  our  results  show  the  prospect 
that  using  advanced  platforms  developed  in  this  project,  analysis  and  prototyping  of  cer¬ 
tain  plasma  device  components  of  interest  to  the  Air  Force  are  amenable  to  cost  efficient 
numerical  simulation. 
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Chapter  2 

Computational  Model 


Here  we  describe  technical  details  of  our  numerical  algorithms  and  implementation  strate¬ 
gies  used  for  physical  modeling  of  high  energy  density  pleismas  of  interest  in  this  project. 
Two  main  codes  will  be  described,  RIGEL  and  ELPIC-MHD,  which  have  been  developed 
here. 

2.1  Physical  Models 

2.1.1  Gyrokinetic/MHD  Hybrid  Scheme  and  Applicability  Lim¬ 
its 

Fluid  Equations 


p^v  +  p(vV)v=  -  Vp-(V-Ph)x  +  JxB  (1) 

— B  =  VxE,  E  =  (vxB-t/J),  J  =  VxB  (2) 

|^  +  V-(,,v)  =  0  (3) 
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—p  +  V- Vp  =  -7pVv  +  />V||  •  K||V||(p/p)  (4) 

VIA  =  -  Jx  ■  (5) 

Gyrokinetic  equations  for  energetic  particles 

Ph  =  Pil  +  {P-  Pi)bb,  (6) 

P||  =  J  cPvF{mu%  =  j  (7) 


F  =  F(R,  u,//) 

dR/dt  =  u[b  +  (u/fi)b  X  (b  •  V)b  ]  +  (l/0)bx(pVR  -  qE/m)  (8) 

du/dt  =  -  [  b  +  (u/fi)b  X  (b  •  V)b  ]  •  {fiWB  -  qE/m)  (9) 

=  -  e(n.  -ne)  (10) 

Ad 

Efficient  implementation  of  the  fast  Laplace  solvers  involved  here  represented  an  im¬ 
portant  part  of  this  effort. 

The  plasma  can  be  modeled  as  single-component  fluid  with  transport  coefficients  inde¬ 
pendent  of  B  under  the  following  conditions.  Let  us  introduce  the  following  notation: 
uj  is  the  typical  flow  frequency  [d/dt  is  of  the  order  of  a;), 

Te,  Ti  are  the  electron-electron  and  the  ion-ion  collision  times,  respectively, 

Qi  is  the  ion  cyclotron  frequency, 

Uj,  Vj  are  the  electron/ion  thermal  velocities,  respectively, 
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va  is  the  Alfven  velocity, 

Pb  is  the  ion  Larmor  radius, 
is  the  ion  Debye  radius, 

Ae,  A,  are  the  mean  free  paths  of  the  electrons  and  ions  respectively, 
lie  =  ne^KmeCo)  is  the  electron  plasma  frequency. 

The  first  condition  necessary  condition  for  the  MHD  approach  to  hold  implies  that  the 
plasma  is  in  a  quasi-stationary  state  characterized  by  a  frequency  u  which  is  much  lower 
■  than  the  effective  frequency  of  electron-ion  collisions: 


1 

-  >  Te,  Ti 

e,t 

(11) 

(jO 

LO 

In  this  case  the  particle  distribution  function  can  be  considered  to  be  isotropic. 

Another 

set  of  conditions 

u>  C  Hi 

i 

^  Pb 

UJ 

(12) 

L  >  min{Xi,pg), 

p'd 

(13) 

results  in  Ohm’s  law  written  cis  follows 

j  =  <T^E-l-vxB 

-  — Vp.)  , 

eUe  / 

(14) 

where  e  is  the  electron  charge,  rzg  is  the  electron  number  density,  and  p,-  is  the  ion  pressure. 
When  |Vpi/n|  ~  |VTi|  <C  \eE\  {Ti  being  the  ions  temperature),  the  Vp,-  term  may  be 
neglected.  Then,  Ohm’s  and  Ampere’s  laws  take  the  following  form 

^  =  E  +  vxB  (15) 

a 

pj  =  V  X  B  (16) 


If  the  following  condition  holds 


then  the  displacement  current  may 


n^Te 


>  1, 


(17) 


be  omitted  in  the  Maxwell  equation  which  reduces  to 


m 

dt 


=  -V  X  E 


(18) 
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The  standard  MHD  field  equation  is  readily  obtained  from  equations  (15),  (16)  and  (18). 
Let  us  now  look  at  specific  figures.  Some  relevant  plasma  parameters  in  mks  units  are^’®: 


n^eHnh.  Ink  \  e  /  VlO^^ 


n.Z^eHnk  ®  Z^lnkKt)  UO^' 


^  ^  g  gg  ^  10^ 

m;  A 


=  4.19  X  lOM  — 


=  1.02  X  10-4 

^m,/  fit  Z  \  e  J  B 

foT,\'l'‘  _ _ V 


A,  =  =2.4x10-= 

\me/ 

A,  =  =  2.4  X  10-S 


7  I'TeV^"  f  n^\ 

(7)  (l^) 

1  /r.\^  / 


InkKeJ  VIO20 


1  fTi 


ZHnk\e)  VIO20 


=  5.64  X  10^ 


m 


Here  Ze  is  the  electrical  charge  of  ions,  A  is  the  atomic  weight  of  ions,  and  Ink  is  the 
Coulomb  logarithm, 

/  n.  \i/2\ 


InA  7  +  2.3logio  f  / 


For  the  plasma  electron  temperature  Te/e  of  10  eV  and  the  electron  number  density  of 
we  obtain  Ink  6,  Tg  ~  2  x  lO^^^s,  r,-  ~  lO^^/Z^s,  fi,-  ~  lO^ZHs  Ag  ~ 
4  X  10“®  m.  A,-  «  4  X  10~^  fZ^*^  m,  Vj  ~  10®  m/s,  Uy  «  3  x  10^  m/s,  ~  2  x  10^  H  m/s, 
fa  3  X  10“^/(Z5)m,  Hg  fa  6  X  lO^^s'h  Substituting  another  typical  set  of  values 
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of  Te/e  =  lOOeV  and  yields  InX  «  8,  Te  «  4  x  r.  w  2  x  IQ-^/Z^s, 

Hi  «  10^  ZB  s~^,  «  10^  m/s,  p'g  «  10~^ {{ZB)  m,  He  ~  2  x  10^“*  s“*.  In  both  cases  all 

the  conditions  (11)  -  (13)  and  (17)  hold  with  confidence. 

2.1.2  The  MHD  approximation 

The  Eulerian  form  of  the  MHD  equations  implemented  in  our  numerical  scheme  is  as 
follows: 

Hydrodynamic  equations 

+  (v  •  V)v  =  --VP  +  i/VV  +  ii/V(V  •  v)  +  -j  X  B  (31) 

ut  p  op 

|^  +  V(,,v)  =  0  (32) 

Magnetohydrodynamic  equations 

^  =  —W^B  T  V  X  (v  X  B)  (33) 

ot  per 

/xj  =  V  X  B  (34) 

V  •  B  =  0  (35) 

where  p  ■=  Air  x  10~^  kg  mjC^  is  magnetic  permeability  of  vacuum,  a  is  electrical  conduc¬ 
tivity  of  plasma. 


2.1.3  Lagrangian  formulation 

The  Lagrangian  equations  implemented  in  our  codes  result  from  a  change  of  variables  in 
the  Eulerian  equations  described  above  into  the  Lagrangian  coordinates.  For  example,  in 
two-dimensional  cylindrical  geometry,  the  transformation  of  the  derivatives 

—I  —I  -I 

(where  quantities  to  the  lower  right  of  the  vertical  bars  are  those  that  are  held  fixed), 
appearing  in  the  Eulerian  equations,  are  obtained  using  the  Jacobian, 
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^  d{r,z,t)  ^  d{r,  z) 
d{R,Z,i')  d{R,Z) 

where  we  have  used  t'  =  t.  The  partial  derivatives  in  the  Eulerian  equations  can  be  written 
as 

d  .  d{  ,z,t)  _  a(H,V,t')  _  1  9{  ,z) 

-  d{r,z,t')  "  “  Jd{R.Z) 

Similarly, 


A|  ^  =  1  9{  ,r) 

dz'"'"  Jd{R,Z) 

The  time  derivative  becomes 


d.  1  5(  ,r,z)  _  1  d{r,z)  d  d(r,z)  d  d{r,z)  d 

m'"’'"  ~  ld{t,R,Z)  ~  J^d{R,Z)dt  d{t,Z)dR  d{t,R)dZ 

Upon  making  the  Lagrangian  identification  dRfdt  =  u  and  dZjdt  =  u,  we  obtain 


dtr,z 


+  UV= 


which  expresses  the  relationship  between  the  Eulerian  and  Lagrangian  time  derivatives 
with  U  =  {u,  v). 


2.1.4  Boundary  conditions  in  the  current-sheet  approximation 

We  have  developed  a  set  of  fast  routines  for  computation  of  flow  regions  where  the  sheet 
current  approximation  holds.  This  speedup  results  from  the  following  physical  considera¬ 
tions.  Since  the  specific  electrical  conductivity  of  plasma  with  Tg  ^  10^  eV  is  of  order  of 
10^(11  m)“^,  the  typical  value  of  of  magnetic  diffusivity  A  =  {ficr)~^  is  of  order  of  Im^fs 
which  implies  that  in  many  physical  setups  of  interest  here  the  magnetic  Reynolds  number 
Re^  =  VL/\  >>  1.  Indeed,  for  plasma  flow  with  characteristic  velocity  of  lO^m/s  and 
characteristic  length  L  of  1  m,  Rcm  is  of  order  of  10®.  This  means  that  only  within  lim¬ 
ited  regions  where  B  changes  significantly  over  a  very  short  scale  8  the  gradients  are  high 
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enough  that  the  magnetic  diffusion  may  not  be  safely  neglected.  In  the  rest  of  the  flow, 
magnetic  field  lines  are  frozen  in  the  plasma  in  a  sense  that  all  fluid  particles  remain  on 
the  field  lines  where  they  initially  were. 

Balancing  the  diffusion  time  5^/A  and  the  convection  time  LjV  results  in  the  estimate 
6  oc  L  *  which  means  that  Sj L  is  of  order  of  10“^  or  less.  In  this  case,  plasma 

is  magnetically  confined  by  current  in  sheet  layer  of  the  width  8  which  appears  on  the 
plasma  boundary.  This  magnetically  confined  plasma  has  a  free  surface  with  no  condition 
on  the  velocity  component  tangential  to  the  boundary  with  vacuum.  The  normal  velocity 
component  Vn  is  equal  to  the  normal  velocity,  if  any,  of  the  interface  itself.  The  condition  on 
the  nor mal-to-the- boundary  component  of  the  magnetic  field,  Bni  requires  that  the  change 
of  Bn  across  the  interface  must  be  zero.  On  the  other  hand,  if  the  plasma  is  considered 
as  perfectly  conducting,  there  may  be  discontinuities  in  tangential  components  of  B.  The 
jump  AB  on  the  interface  implies  a  surface  sheet- current.  By  Stokes’  theorem  applied  to 
equation  (34),  the  sheet-current  density  per  unit  distance  measured  tangentially  normal  to 
the  current  is 

J  =  -n  X  AB,  (36) 

IJ' 

n  being  a  unit  vector  normal  to  the  interface.  Of  course  if  plasma  is  considered  as  finitely 
conducting,  the  tangential  field  components  must  also  be  continuous  across  the  interface. 
The  extra  condition  of  a  zero  AB  matches  the  rise  in  the  order  of  the  field  equation  (33) 
when  the  magnetic  diffusion  term  is  included.  The  magnetic  boundary  layer  of  a  finite 
thickness  appears  in  this  formulation.  It  can  be  viewed  upon  as  a  dilated  version  of  a 
sheet-current  made  diffuse  by  the  finite  conductivity. 

This  magnetic  field  jump  condition  implies  certain  force  balance  at  the  plasma  interface, 
namely,  a  discontinuity  of  pressure  across  the  current  sheet  within  which  the  normal  j  x  B 
force  acts.  Since  the  plasma  is  confined  by  its  interface,  the  fluid  acceleration  is  finite  and 
the  inertia  term  does  not  appear  in  the  jump  conditions.  The  magnetic  force  per  unit 
volume  consists  of  two  terms; 

j  X  B  =  -(V  X  B)  X  B  =  -  ((B  •  V)B  -  VBV2)  (37) 

fi  fj, 
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Inside  the  sheet-current,  the  first  contribution  (B  •  V)B//i  is  estimated  as  where  R 

is  the  curvature  radius  of  the  interface.  The  second  term  VB^/2/i  is  of  order  of  j fxS  and 
obviously  dominates.  Therefore,  the  force  balance  condition  at  the  boundary  is 

A{p  +  By2fi)  =  0,  (38) 

where  only  the  tangential  components  of  B  contribute. 

2.1.5  Plasma  properties 

We  have  developed  an  extensive  database  of  material  properties  such  as  source  terms  for 
radiation  transport  and  equations  of  state  (EOS).  The  central  issue  in  modeling  of  radiation 
transport  is  implementation  of  libraries  of  different  models  describing  the  emission  and  ab¬ 
sorption  of  radiation  which  are  required  depending  on  the  simulation  under  consideration. 
Local  thermodynamic  equilibrium  (LTE)  models  are  adequate  for  modeling  radiation  trans¬ 
port  in  the  high-density  colder  regions  plasmas  while  a  full-LTE  rate-equilibrium  model  is 
used  to  accurately  model  radiative  processes  in  and  around  hot,  lower-density,  high-atomic- 
number  plasmas.  In  many  simulations,  the  coupling  of  the  radiation  to  the  fluid  is  treated 
explicitly.  The  energy  removed  from  or  deposited  to  the  electrons  by  each  energy  group  is 
summed  over  all  energy  bins  transported  for  each  computational  zone.  The  source  term  is 
then  explicitly  added  to  the  electron  thermal  energy  equation.  In  some  cases,  however,  the 
radiation  and  fluid  are  closely  coupled  and  must  be  solved  as  such. 

The  analytic  Thomas-Fermi/SESAME  EOS  was  derived  for  use  in  instability  simula¬ 
tions  here.  Tabular  EOS’s  such  as  SESAME,  while  more  accurate  than  the  analytic  EOS, 
have  been  found  to  result  in  excessive  numerical  noise  during  the  simulation,  to  the  extent 
of  reaching  high  enough  levels  to  render  the  entire  simulation  invalid.  For  the  non-ideal 
gas  EOS  options,  we  have  implemented  a  tabular  look-up  of  the  average  z  and  average  . 
These  tables  are  based  upon  information  obtained  from  the  2000-group  LTE  opacity  tables 
used  in  the  radiation  transport  modules  of  our  codes. 
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2.2  Numerical  Implementation  Aspects 

2.2.1  Treatment  of  the  Artificial  Viscosity  Term  (AV) 

The  following  general  requirements  were  targeted  for  incorporation  of  the  artificial  viscosity 
in  our  code  modules 

(1)  Differential  equations  (  with  a  dissipative  mechanism  included  )  should  apply  ev¬ 
erywhere  in  the  interior  of  flow,  with  no  internal  boundary  condition  required; 

(2)  the  basic  conservation  laws  expressed  by  the  Rankine-Hugoniot  conditions  should 
be  obtained; 

(3)  shocks  should  manifest  themselves  as  approximate  discontinuities  in  corresponding 
fields; 

(4)  the  thickness  of  a  smeared  shock  should  be  independent  of  the  shock  strength  and 
the  material  through  which  the  shock  travels; 

(5)  the  artificial  viscosity  should  be  independent  of  homogeneous  expansions  or  contrac¬ 
tions  of  the  medium  (  that  is,  in  general,  thermodynamically  reversible  processes  should 
not  give  rise  to  shock  formation); 

(6)  the  velocity  component,  parallel  to  the  shock  front  in  a  medium,  should  be  contin¬ 
uous; 

(7)  angular  momentum  should  be  conserved;  and 

(8)  there  should  be  no  artificial  viscosity  for  a  velocity  field  under  rigid  rotation. 

A  number  of  different  expressions  have  been  developed  and  implemented  into  our  mul¬ 
tidimensional  simulation  codes.  One  of  the  best  user  options  has  the  important  feature  of 
reducing  to  the  Richtmyer  and  von  Neumann  form  for  one-dimensional  flows  for  both  two- 
and  three-dimensional  implementations: 

AV  =  alp(8v)^,Sv  <  O',  (39) 

AV  =  0,Sv>0. 

Here  p  is  the  fluid  density,  8v  is  the  change  in  velocity  across  a  zone,  and  is  a  constant 
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related  to  the  number  of  computational  mesh  points  over  which  the  shock  is  smeared.  This 
form  of  AV  satisfies  most  of  the  desired  properties.  In  particular,  Condition  2  is  satisfied, 
provided  the  dimensions  over  which  the  shock  is  smeared  are  small  compared  to  other 
dimensions. 

2.2.2  Numerical  Treatment  of  Boundary  Conditions 

The  treatment  of  boundary  conditions  for  the  hydrodynamics  in  our  Lagrangian  code  mod¬ 
ules  is  accomplished  without  any  special  treatment  of  the  momentum  equation  at  the 
boundary.  Nonrigid  boundaries  are  treated  by  assuming  the  pressure  is  zero  (or  prescribed) 
in  the  region  outside  the  computational  mesh. 

The  proper  boundary  condition  for  thermal-energy  diffusion  is  the  zero  gradient  (“no 
flux”)  boundary  condition  at  a  vacuum-fluid  interface.  To  accommodate  this  no-flux  con¬ 
dition,  it  is  necessary  to  extend  the  computational  grid  when  computing  temperature  gra¬ 
dients  near  a  vacuum  interface  (“ghost  cell”).  The  ghost  cells  are  used  only  to  enforce  the 
boundary  condition  and  do  not  affect  material  zones  in  the  problem.  At  slip  surfaces,  the 
zero-temperature  gradient  condition  is  also  implemented. 

2.2.3  Thermal  flux  limitation 

It  hcis  been  shown  by  a  number  of  authors  that  the  Spitzer-Harm  electron  thermal  conduc¬ 
tivity  used  in  most  high-energy-density  simulations  breaks  down  under  two  conditions. 
The  first  is  at  initial  low  “room”  temperature  and  normal  density  where  the  Spitzer-Harm 
conductivity  will  result  in  low  values.  The  second  is  in  regions  where  the  electron  mean 
free  path  is  large  compared  to  the  thermal-gradient  scale  length.  In  the  latter  case  Fouriers 
law  for  thermal  conduction  is  invalid.  In  this  work,  we  followed  the  existing  practice  in 
high-energy-density  plasma  simulation  research  to  limit  the  thermal  flux  to  a  maximum 
value  corresponding  to  the  energy  transported  by  free-streaming  particle  flow^^.  There  are 
two  methods  commonly  used  for  this  limiting  procedure.  In  both  cases  an  effective  free- 
streaming  conductivity  is  found  by  dividing  the  free-streaming  flux  by  the  temperature 
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gradient.  In  the  first  case  the  conductivity  used  in  the  conduction  equations  is  set  equal 
to  the  minimum  value  between  the  Spitzer-Harm  thermal  conductivity  and  the  effective 
free-streaming  conductivity.  The  second  method  is  to  use  a  harmonic  mean  between  the 
two  conductivities,  thereby  allowing  for  a  smooth  transition  between  the  thermal  conduc¬ 
tion  given  by  Spitzer-Harm  and  the  free-streaming  value.  The  application  of  a  flux  limiter, 
while  physically  motivated,  is  ad  hoc.  The  use  of  flux  limiters  results  in  approximately 
physically  correct  behavior. 

2.2.4  Grid  distortion  control  and  rezoning  techniques 

A  number  of  algorithms  have  been  implemented  in  our  code  modules  in  order  to  mini¬ 
mize  the  development  of  nonphysical  grid  distortions.  For  both  the  two-and  three-  dimen¬ 
sions  options,  the  tensor  artificial  viscosity  of  stabilizing  grids  similar  to  that  developed 
by  Wilkins  ^  has  been  implemented.  In  the  two-dimensional  modules,  an  additional  grid- 
distortion-control  option  has  been  implemented  based  on  a  rotational  artificial  viscosity 
originally  developed  for  the  tensor  code 

It  should  be  mentioned  that  rezoning  of  one-,  two-and  three-dimensional  Lagrangian 
hydrodynamics  programs  is  becoming  an  important  area  of  research  in  code  development 
in  the  high-energy-density  plasma  simulation  field.  Rezoning  represents  the  Eulerian  step 
associated  with  the  arbitrary  Lagrangian/Eulerian  implementation  of  the  hydrodynamics 
routines  in  RIGEL.  During  this  phase  of  the  simulation  the  computational  mesh  can  be 
rezoned  (either  in  a  fully  Eulerian  manner  or  by  some  other  procedure)  and  all  advective  flux 
calculations  are  performed.  When  flow  distortions  occur  that  result  in  unacceptably  small 
time  steps  or  poor  resolution  in  regions  of  interest,  some  corrective  action  is  required.  In 
our  code  modules,  we  use  several  rezoning  algorithms,  including  van  Leer’s  methods 
Donor-cell  differencing  has  been  implemented  for  the  one-,  two-,  and  three-  dimensional 
options. 


12 


2.2.5  Technical  aspects  of  solver  implementation 

Here  we  illustrate  some  issues  and  methods  used  in  discretization,  solvers,  and  parallel 
implementation.  Consider  for  example  the  energy  diffusion  in  two  dimensions  on  a  distorted 
mesh  which  also  includes  the  effects  of  heat  flow  from  the  “corner”  zones.  Corner  zones 
are  defined  as  those  surrounding  the  zone  of  interest  but  not  directly  adjacent  (sharing  a 
common  interface)  to  that  zone.  The  method  used  in  this  work  builds  upon  our  algorithm 
previously  developed  for  the  ORCHID  code  and  consists  of  a  combination  of  both  finite 
difference  and  methods  more  closely  related  to  finite  elements^  Let  us  designate  the 
temperature  at  the  intersection  point  of  branch  line  j,  zone  i,  as  Let  also  Tc^.  designate 
the  temperature  defined  at  the  area  centroid  of  zone  k  so  that,  say,  the  temperature  at 
vertex  point  1  is  obtained  by  averaging  of  centroid  temperatures 

Tci  +  Tc2  +  Tcj  +  Tc^ 

-'vi  -  ^ 

Approximation  of  the  temperature  gradient  involves  the  following  length  scales:  A  is  the 
length  of  the  panel  from  vertex  point  1  to  point  2,  Ai  is  the  length  along  the  panel  from 
vertex  point  1  to  the  intersection  point  of  the  normal  to  the  panel  from  the  area  centroid 
of  zone  1;  A2  is  the  length  along  the  panel  from  vertex  point  1  to  the  intersection  point  of 
the  normal  to  the  panel  from  the  area  centroid  of  n  zone  2;  Axj  2  are  the  normal  distances 
from  the  branch  line  intersection  points  in  zones  1,2  to  the  panel,  respectively. 

The  form  of  the  temperature  gradients  depends  upon  which  branch  lines  the  panel 
perpendicular  bisector  intersects.  The  expressions  for  temperature  gradient  for  the  four 
possible  cases  follow: 


(a)  For  Ai  <  A/2  and  A2  <  A/2, 


3:i(A/2) 

A-Ai  ’ 


Aa:2 


^2(A/2)  _  (A/2  -  A0r„2  +  (V2)r,, 

A  -  A2  ’  A  -  Ai 
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VT  = 


Axl  +  Aa;} 


(6)  For  Ai  >  A/2  and  Aj  >  A/2, 


A  _2  _  ^1(^/2)  A  „2  _  ^2(A/2)  ^2  _  (^1  “  ^I‘^)Tv2  +  (A/2)Tt 


Al 


/7-t2  _ 

■’  ^61  - 


Cl 


Ai 


(A,  -  A/2)r,.  +  (A/2)r„ 
'»>  -  1. 


„  W.  -  r^) 

■  =  A.1  +  Axj- 


The  two  other  possible  cases  are 


(c)  Aj  >  A/2  and  A2  A/2, 


(d)  Ai  <  A/2  and  A2  >  A/2. 

Similar  derivation  using  the  correct  form  of  intersection  point  temperature  and  the 
normal  from  the  intersection  point  to  the  panel  yields  the  expressions  for  the  temperature 
gradient  for  these  final  two  cases: 


W  VT 


{Ti  -  n,) 

Ax\  +  Axf’ 


^  {Tl  -  Tl  ^ 

Arcf  +  Axi 

Although  the  temperature  gradient  still  retains  the  simple  form  AT / Ax  it  is  nonlocal 
and  accounts  for  both  grid  distortion  and  the  influence  of  surrounding  ^^corner”  zones.  This 
procedure  results  in  a  sparse  matrix  which  is  inverted  using  the  successive  over-relaxation 
(SOR)  variation  of  the  Gauss-Seidel  method.  Consider  a  set  of  m  linear  algebraic  equations, 
written  in  matrix  form  as 
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^n+1  ^  ^  +  UX^  +  h  -  X”), 

where  the  superscript  denotes  the  iteration  number.  The  change  in  x  between  two  iterations 
(^n+i  _  ig  referred  to  as  the  correction  to  x  or  the  displacement  vector.  The  relaxation 
parameter  lo  is  usually  a  real  positive  number  between  0  and  2.0.  For  uj  less  than  1  the 
method  termed  is  “under- relaxed” .  Note  that  w  =  1  reduces  the  SOR  iteration  to  Gauss- 
Seidel  iterates. 

Both  the  Gauss-Seidel  and  SOR  methods  replace  x"  by  .t"'*''  as  soon  as  becomes 
available.  SOR  is  performed  using  a  relative  error  convergence  criterion;  x"'*’^  —  x")/x"  <  e, 
where  e  is  chosen  typically  in  the  range  10“®  to  10“®  depending  on  the  problem. 

On  standard  single  to  “few”  central-processor  implementations  we  have  implemented 
SOR  based  on  the  Gauss-Seidel  iteration  procedure  whereupon  determining  and  updating 
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iterate  it  is  used  in  succeeding  calculations  of  other  iterates.  For  a  “many”-  processor 
implementation  there  is  an  option  for  a  Jacobi-type  iteration  where  the  iterates  are  updated 
only  at  the  end  of  the  overall  system-iteration  step  (sweep).  This  solution  procedure,  while 
easier  to  implement  on  a  many-processor  computer  system,  represents  an  area  of  further 
research  and  development.  The  procedures  of  determination  of  the  gradient  on  a  distorted 
mesh  in  two  dimensions  and  three  dimensions  are  similar.  In  three  dimensions,  instead  of 
intersections  of  the  panel  normal  with  linear  interpolation  segments  between  nodal  points 
and  zone  centroids,  intersections  with  interface  normal  and  planes  is  used,  resulting  in  a 
27  diagonal  system. 

The  current  version  of  ELPIC-MHD  includes  parallelization  options  which  force  the 
most  time-consuming  internal  loops  of  the  code  run  in  parallel  on  many  processors.  The 
advantage  of  using  a  parallel  code  becomes  apparent  taking  into  account  that  95  percent  or 
more  of  the  entire  CPU  time  in  computation  of  the  magnetic  field  is  consumed  by  the  inner 
nested  loops  which  implements  Gauss-Seidel  algorithm  with  overrelaxation  (SOR).  In  this 
loop,  the  central  processor  recalculates  full  arrays  of  values  of  the  magnetic  field  component 
along  with  the  current  local  iteration  residual  matrix  R  as  well  as  the  average  residual  for 
the  entire  matrix.  At  each  time  step,  up  to  several  thousands  of  system-iteration  sweeps 
may  be  required  for  each  of  the  magnetic  field  components.  On  a  single  processor  this 
procedure  routinely  takes  at  least  70-80%  of  the  total  CPU  time  on  large  grids  of  practical 
interest. 

In  the  hydrodynamical  part  of  the  code,  transfer  of  Lagrangian  particles  takes  60  to 
75%  of  the  rest  of  the  CPU  time.  The  core  part  of  this  procedure  is  coded  as  a  loop  of 
about  600  executable  statements  which  consumes  almost  100%  of  the  total  particle-pushing 
CPU  time.  This  loop  is  executed  once  for  each  one  of  the  Lagrangian  particles. 

Together,  these  two  loops  take  90  to  over  95%  of  the  total  CPU  time.  Therefore,  the 
maximum  acceleration  by  as  much  as  15-20  can  be  achieved  by  optimal  parallelization 
of  only  these  two  loops  on  20-30  processors.  Parallelization  of  the  Lagrangian  stage  is 
straightforward.  For  a  “few”  processors,  the  most  efficient  approach  is  to  use  the  par- 
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allelization  directives  at  the  source  level.  Source-independent  parallelization  algorithms 
have  been  found  to  split  the  loops  in  such  a  way  that  particle  pushing  affects  distant  ele¬ 
ments  of  the  corresponding  Eulerian  grids,  automatically  providing  a  high  efficiency  of  such 
parallelization. 

The  bottleneck  of  speedup  due  to  optimization  of  the  Eulerian  stage  using  source-level 
directives  is  efficient  treatment  of  second  order  differential  operators.  The  stencil  here 
involves  several  elements,  some  of  which  result  from  calculation  using  the  same  formula 
obtained  earlier  in  the  nested  loop.  This  dependency  is  required  to  take  advantage  of 
overrelaxation  formula  which  helps  to  dramatically  reduce  the  total  number  of  iterations. 
Full  explicit  parallelization  of  this  stage,  performed  using  source-level  directives  which 
divide  loops  into  several  parts  and  force  the  processors  to  execute  those  parts  in  parallel,  is 
not  efficient.  The  reason  is  that  the  processors  have  to  wait  until  the  values  of  the  function 
along  the  adjacent  subdomains  are  calculated  by  the  neighboring  processors. 

An  alternative  parallelization  strategy  provides  almost  linear  scaling  of  performance 
with  the  number  of  processors.  The  idea  here  is  as  follows.  Convergence  criteria  outside  the 
stencil  at  given  point  in  the  computational  domain  are  always  known  while  the  operations 
at  this  given  point  are  performed.  At  the  same  time,  at  least  dozens  of  sweeps  are  performed 
over  the  same  array  within  a  time  step.  The  optimized  algorithm  scans  the  array  line-by- 
line  starting  from  the  line  next  to  the  domain  boundary.  After  scanning  of  three  lines  is 
completed  by  the  first  processor,  the  same  loop  starts  on  the  second  processor  which  begins 
scanning  again  from  the  first  line.  The  array  elements  read  and  written  by  the  second 
processor  axe  fully  independent  of  the  elements  read  and  written  by  the  first  one,  since  the 
calculated  elements  are  separated  by  at  least  two  lines  of  the  array.  Computational  stencils 
involved  in  reading  and  writing  operations  do  not  overlap.  Therefore,  two  consequent 
internal  iterations  are  performed  in  parallel  by  two  individual  processors  with  virtually  no 
overhead.  When  the  second  processor  completes  calculation  of  the  first  three  lines,  the 
third  processor  takes  up  the  job  beginning  with  the  first  line  and  so  on.  When  the  first 
processor  completes  execution  of  the  last  line  of  the  matrix,  it  checks  whether  the  last 
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processor  has  completed  execution  of  the  first  three  lines  and  if  this  condition  is  met  it 
starts  scanning  the  first  line  again.  All  the  other  processors  follow  the  same  algorithm. 
After  a  preset  number  of  full  sweeps,  one  of  the  processors  (in  our  implementation,  the  last 
one)  checks  the  overall  convergence  criteria  for  the  entire  matrix  using  the  lower  estimate 
of  the  residual  norm  which  it  accumulated  during  the  last  sweep  (the  real  residual  norm 
can  only  be  less).  If  the  norm  does  not  exceed  the  known  (preset  or  calculated)  value,  the 
loop  execution  is  completed. 

In  this  way,  using  up  to  N/3  -  1  processors  where  N  is  the  total  number  of  lines  result  in 
minimal  degradation  of  relative  performance.  This  cyclic  algorithm  provides  almost  linear 
scaling  if  the  required  total  number  of  loops  is  much  greater  than  the  number  of  available 
processors.  Using  the  limited  depth  multigrid  method  improves  performance  even  further. 
If  the  number  of  multigrid  levels  is  limited  by  two,  the  algorithm  is  identical  to  the  one 
described  above,  with  the  only  difference  that  the  number  of  processors  which  can  be 
effectively  used  in  parallel  may  not  exceed  M/3-1,  where  M  is  the  number  of  lines  in  the 
coarser  subgrid. 

In  simulations  involving  “very  many”  processors,  the  arrays  containing  particle  data 
have  to  be  sorted  out  so  that  particles  from  consequent  elements  of  these  arrays  always  be¬ 
long  to  the  cells  of  Eulerian  grid  separated  by  at  least  two  other  cells.  For  three-dimensional 
grids  this  critical  number  of  processors  is  about  (1/3)^  w  1/30  of  the  total  number  of  Eule¬ 
rian  cells  which  ranges  from  tens  of  thousand  to  millions  depending  upon  the  application. 
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Chapter  3 


Numerical  Modeling  of  Plasma 
Devices 

Our  computational  strategies  described  in  the  previous  chapter  can  be  applied  to  a  wide 
variety  of  practically  important  issues  such  as 

-  advanced  weapons  systems; 

-  pulsed  plasma  radiation  sources; 

-  advanced  propulsion  systems,  including  prospective  space-application  thrusters  based 
on  burning  of  D  —  He^,  e.g.  nuclear  thermal  propulsion  thrusters; 

-  microwave  generation; 

-  inertial  confinement  fusion. 

In  this  chapter  we  review  several  prototype  simulations  of  high  energy  density  plasma 
devices  of  interest  for  the  AirForce,  emphasizing  plasma  focusing  and  energy  conversion. 
These  test  cases  are  representative  of  a  variety  of  parametric  studies  that  we  performed  in 
the  course  of  this  project,  including  variable  device  geometries,  physical/material  proper¬ 
ties,  boundary/initial  conditions,  different  code  modeling  modules  options  used,  etc. 

The  first  test  case  presented  here  involves  a  prototype  configuration  of  a  plasma  beam 
compression  device  used  for  pulse  magnetic  field  confinement  of  injected  plasma  beams.  In 
this  device,  compression  is  due  to  the  configuration  of  magnetic  field  in  the  vacuum  chamber 
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of  nearly- conical  shape.  The  device  responsible  for  generation  of  short  pulse  magnetic  field 
is  located  outside  our  computational  domain.  The  initial  and  boundary  conditions  in  the 
runs  illustrated  in  Figs  1-18  correspond  to  relatively  cool  fully  ionized  homogeneous  plasma 
with  a  uniform  coefficient  of  ionization  injected  through  the  inlet  located  on  the  right  side 
of  each  plot.  The  beam  is  formed  and  pre-conditioned  during  homogenization,  acceleration, 
and  preliminary  compression  stages.  The  injection  speed  used  in  presented  runs  is  50  km/s 
(the  entire  range  of  injection  speeds  from  10  to  150  km/s  has  been  studied).  The  inlet 
plasma  density  in  the  beam  drops  sharply  from  a  typical  value  of  O.OSif^r/m^  in  the  core 
to  very  small  values  outside  of  a  tube  of  about  75%  of  the  overall  inlet  radius.  Most  of  the 
magnetic  flux  injected  through  the  inlet  is  concentrated  in  the  area  of  low  plasma  density 
restricted  on  the  outside  by  the  inner  surface  of  the  device  wall.  The  magnetic  flux  is  set 
constant  over  duration  of  plasma  injection.  The  integral  magnetic  fluxes  through  the  inlet 
and  through  the  outlet  (located  on  the  left  side  of  our  plots)  are  equal. 

Since  the  injected  beam  is  confined  by  magnetic  field,  dense  plasma  is  isolated  from  the 
outer  wall.  The  possibility  to  set  the  tangential  component  of  particle  velocity  at  the  wall 
to  zero  is  embedded  in  our  code  modules  in  order  to  properly  model  wall  boundary  layer 
effects  important  for  device  efficiency  and  durability  analysis. 

This  magnetically  confined  plasma  beam  propagates  through  the  chamber  with  smoothly 
diminishing  cross-section  and  incrccising  density.  Our  setup  allows  for  a  parametric  study 
of  the  behavior  of  this  magnetically  confined  and  compressed  plasma,  thus  prototyping 
a  practically  important  experimental  study  aimed  at  achieving  maximum  density  in  the 
plasma  jet  near  the  outlet  accompanied  by  maximum  acceleration/inhibited  deceleration 
and  therefore  avoiding  inasmuch  as  possible  conversion  of  the  kinetic  energy  into  heat.  This 
optimization  is  important  in  several  applications  which  involve  an  outcoming  jet  hitting  a 
target  in  order  to  generate  x-rays  pulses  of  desired  characteristics  with  maximum  intensity. 

In  Figs.  1-7  we  show  a  complete  time  history  of  a  run  in  which  the  magnetic  field 
configuration  is  near  the  optimum  from  the  point  of  view  of  above  considerations.  In  Fig. 
1  we  plot  the  initial  setup  with  the  magnetic  field  lines  plotted  in  green  color. 
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Initially  the  chamber  is  nearly  empty.  A  highly  conducting  plasma  beam  is  injected 
through  the  inlet  on  the  right  side  of  the  picture.  The  specific  internal  energy  of  different 
plasma  regions  is  shown  in  terms  of  a  color  palette  imitating  those  of  a  heated  body. 
Dark  grey  color  stays  for  cold  material,  whereas  shades  of  red  transforming  into  yellow, 
then  green,  blue,  and  at  last  white  indicate  increasing  energy  and  temperature.  The  color 
scale  on  top  of  Fig.  1  shows  the  color-energy  correspondence.  For  convenience,  the  specific 
internal  energy  is  expressed  in  degrees  Kelvin. 

In  the  lower  part  of  Fig.  1,  we  plot  the  very  beginning  of  plasma  compression.  The 
effective  magnetic  field  pressure  applied  to  the  plasma  surface  generates  a  shock  wave.  At 
the  same  time,  the  front  axial  part  of  the  beam  does  not  experience  any  significant  magnetic 
pressure  and  tends  to  slightly  expand  and  cool  down. 

The  compression  history  is  shown  in  Fig.  2.  The  compression  shock  at  the  surface  which 
is  initially  strong  is  partly  relieved  by  a  slight  expansion  of  the  frontal  surface  of  the  beam. 
This  expansion  is  faster  in  the  core  area  where  the  effective  magnetic  pressure  is  minimal, 
producing  a  clearly  visible  distortion  of  the  frontal  surface  of  the  plasma  beam  (which  was 
initially  set  flat). 

In  Fig.  3  we  show  the  most  interesting  stage  of  the  compression  process,  the  focusing. 
In  the  bottom  part  of  Fig.  3,  the  very  moment  of  focusing  is  shown.  Shortly  before  the 
focusing  event,  the  surface  shock  reconnects  along  the  entire  plasma  surface.  The  bright 
spot  shows  the  exact  position  of  the  focal  point.  Until  the  very  moment  of  focusing  the 
plasma  velocity  slowly  decreases,  but  right  after  the  focusing  the  plasma  beam  tip  on  the 
left  side  of  the  focal  point  starts  to  accelerate. 

In  Fig.  4  we  plot  developing  hot  area  in  the  beam  core  produced  by  collapse  of  the 
imploding  shocks  on  the  axis,  resembling  that  in  Munro  jets.  The  velocity  of  forming  axial 
hot  and  dense  jet  is  always  directed  against  the  mean  flow  velocity.  Therefore,  the  kinetic 
energy  of  the  initial  plasma  beam  is  partly  converted  into  heat  and  tends  to  slow  down  the 
bulk  flow.  Excessive  core  heating  poses  a  limitation  on  the  maximum  compression  rates 
that  can  be  achieved. 
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In  Fig.  5  one  can  see  formation  of  a  metastable  structure  in  the  focused  plasma  beam 
which  consists  of  a  sequence  of  domains  with  lower  or  greater  thickness.  The  cause  is  clearly 
visible  in  the  top  part  of  Fig.  5:  two  consequent  axial  jets  are  formed,  both  slowing  down 
the  compressed  beam. 

In  Fig.  6,  we  plot  a  developed  axial  jet  that  has  propagated  far  upstream,  progressively 
slowing  down  the  beam.  The  maxima  of  compression  ratio,  temperature,  and  pressure, 
as  well  as  the  velocity  minimum  are  achieved  in  the  area  adjacent  to  the  outlet,  where 
plasma  effectively  expels  the  magnetic  field  from  the  channel.  The  configuration  shown 
in  the  bottom  portion  of  Fig.  5  precedes  a  moment  when  substantial  number  of  particles 
reach  the  wall  blocking  the  way  to  the  magnetic  flux. 

This  blocked  configuration  is  shown  in  Fig.  7.  In  real  life  laboratory  experiment  this 
stage  corresponds  to  magnetic  field  being  completely  expelled  from  the  chamber  so  that  its 
energy  is  almost  instantly  emitted  in  form  of  electromagnetic  waves.  Most  of  the  kinetic 
energy  of  our  initially  cold  incoming  jet  is  immediately  converted  into  shock  heating  of 
the  stagnated  dense  plasma  cloud  near  the  outlet,  in  the  leftmost  part  of  the  chamber 
(bright  white  color  which  depicts  temperatures  in  the  Kev  range).  A  magnified  part  of  the 
stagnation  zone  is  shown  in  the  bottom  portion  of  Fig.  7.  A  very  hot  rarefied  plasma  region 
appears  as  separate  particles  and  propagates  far  upstream  in  the  gap  between  the  jet  and 
the  wall.  These  small  portions  of  plasma  merge  back  with  the  surface  of  the  core  jet  and 
substantially  heat  up  the  thin  surface  layer. 

In  Fig.  8  we  show  a  configuration  corresponding  to  the  same  initial  and  boundary 
condition  of  the  flow,  with  a  magnetic  field  approximately  twice  as  high  as  in  previous  run. 
A  naive  assumption  would  be  that  increased  magnetic  field  should  lead  to  increased  beam 
compression  and  therefore  improved  efficiency. 

However,  data  plotted  in  Fig.  8  indicates  that  such  naive  assumptions  can  be  wrong. 
Indeed,  in  this  run  a  stronger  magnetic  field  initially  makes  the  first  focal  point  closer  to 
the  inlet  and  formes  a  thinner,  denser  compressed  jet.  This  jet  then  accelerates  to  a  slightly 
higher  velocity.  However,  the  rarefied  and  fastest  leftmost  edge  of  the  jet  experiences  rapid 
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deceleration  by  the  magnetic  field.  This  creates  a  shock  at  the  jet  tip,  overheats  it  and 
makes  it  rapidly  expand.  This  in  turn  formes  a  cavity-like  configuration  of  the  magnetic 
field,  leading  to  further  increase  of  plasma  deceleration.  This  positive  feedback  instability 
leads  to  rapid  formation  of  a  dense  hot  stagnated  plasma  cloud.  The  kinetic  energy  is  lost 
into  heat  and  emitted  prematurely  inside  the  chamber.  This  result  strongly  indicates  that 
optimum  plasma  focusing  can  be  achieved  within  nontrivial  parameter  regimes  which  elude 
common  sense  but  require  systematic  numerical  study. 

In  Figs.  9  -  13,  we  plot  a  time  history  of  another  run  characterized  by  near-the-optimum 
magnetic  field  magnitude.  In  these  Figures,  walls  and  magnetic  field  lines  are  not  shown. 
This  simulation  intentionally  starts  from  a  domain  already  filled  with  moving  plasma, 
except  for  a  small  two  cell-wide  gap  between  the  plasma  and  wall.  The  goal  is  to  compu¬ 
tationally  determine  the  minimum  magnitude  of  the  magnetic  field  which  separates  dense 
plasma  from  the  wall  before  it  blocks  the  magnetic  flux.  Fig.  10  shows  the  plasma  config¬ 
uration  at  such  magnetic  field.  This  is  a  nearly  ideal  focusing.  Figs.  11-13  show  the  rest 
of  the  numerical  experiment.  Initially  (Fig.  11)  the  jet  diverges  after  the  focal  point.  At  a 
later  time  (Fig.  12),  the  plasma  touches  the  wall  forming  a  strong  shock.  At  last  (Fig.  13) 
the  stagnation  cloud  blocks  the  channel  near  the  outlet. 

In  Fig.  14  we  show  the  outcome  of  a  run  with  the  same  initial  configuration  but  an 
even  lower  magnitude  of  magnetic  field.  Here,  we  plot  the  walls  but  not  the  magnetic  field 
which  is  concentrated  in  the  gap  next  to  the  wall  and  within  it.  In  this  experiment,  no  focal 
point  appears  at  all.  Shortly  after  the  start,  plasma  contacts  the  wall  forming  a  shock  that 
separates  far  from  the  outlet.  The  shock  collapses  onto  the  axis  and  forms  a  stagnation  core 
which  then  absorbs  and  converts  all  the  jet  kinetic  energy.  Fig.  14  illustrates  the  instant 
of  stagnation  core  formation. 

We  have  seen  that  the  parameter  domain  corresponding  to  absence  of  premature  cham¬ 
ber  blocking  can  be  multiconnected.  In  Figs.  15  -  18  we  plot  a  time  history  of  a  run  with 
the  magnetic  field  amplitude  just  slightly  stronger  than  the  optimum  (c.f.  with  the  series 
shown  in  Fig.  8  which  corresponds  to  only  slightly  lower  magnetic  field).  The  configuration 
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shortly  after  the  start  is  shown  in  Fig.  15.  As  one  can  see  from  Figs.  16-18,  a  dense  jet 
with  a  sequence  of  focal  points  formes  and  exists  for  an  extended  time.  Comparing  Figs.  17 
and  18,  one  can  see  that  the  first  focal  point  retains  its  position  during  the  run  while  the 
second  focal  point  slowly  moves  upstream,  away  from  the  outlet. 

Special  attention  was  given  throughout  the  project  to  the  problem  of  compact  toroid 
acceleration  by  magnetic  field.  The  compact  toroid  is  a  nearly  force-free  plasma  configu¬ 
ration  whose  acceleration  is  provided  by  a  pulse  of  toroidal  magnetic  flux  injected  radially 
inwards.  This  pulse  also  forces  reconnection  of  magnetic  field  lines  to  form  a  closed  toroidal 
plasma-field  object  with  extended  lifetime. 

A  set  of  specialized  subroutines  has  been  specifically  developed  in  order  to  fast  set 
up  nearly  force-free  toroidal  configuration  in  complex  device  geometries.  In  Figs.  19  - 
22  we  plot  the  acceleration  stage  of  compact  toroid  in  a  prototype  device  similar  to  the 
formation  and  compression  regions  of  the  MARAUDER  facility.  In  Fig.  19  we  plot  the  stable 
flow/magnetic  field  configuration  as  created  by  our  specialized  setup  routine.  The  formation 
region  of  the  prototype  device  has  constant  width  and  the  narrowing  zone  corresponds  to 
the  compression  region.  A  nearly  round  cross-section  of  the  compact  toroid  is  located  next 
to  one  of  the  inlets  designed  for  injection  of  a  pulse  acceleration  flux.  In  a  real  device, 
this  corresponds  to  the  zone  near  the  the  muzzle  of  the  formation  gun.  A  blowup  of  the 
outside  cross  section  of  the  magnetic  field  configuration  calculated  by  fast  solver  is  plotted 
in  Fig.  20.  Inside  the  toroid,  equilibrium  distribution  based  upon  balance  of  magnetic  and 
hydro  pressure  is  computed.  Temperature,  density  and  hydro  pressure  smoothly  increase 
towards  the  axis  of  the  toroid,  compensated  by  the  magnetic  effective  pressure  gradient. 
The  total  toroid  mass  is  0.72  x  10“^  grams.  The  total  initial  number  of  computational 
Lagrangian  particles  in  the  toroid  is  205,000.  Since  the  particles  breed  in  the  areas  of  high 
stretching,  their  number  increases  during  the  run. 

The  configuration  at  a  later  time  is  shown  in  Fig.  21.  The  magnetic  field  around  the 
toroid  has  significantly  changed.  Notice  high  field  gradients  near  the  rear  surface  of  the 
toroid.  The  magnetic  force  acting  at  the  rear  surface  is  now  much  above  that  at  the  front 
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surface,  providing  for  the  acceleration  force. 

A  representative  later  stage  of  the  toroid  evolution  dynamics  is  plotted  in  Fig.  22.  The 
magnetic  field  is  not  confined  within  the  toroid  —  field  lines  leave  the  plasma  region  and  then 
return  back  in  the  toroid  body.  The  development  of  small  scale  surface  structures  is  due  to 
the  weak  instabilities.  Significant  reconfiguration  of  the  internal  field  in  the  compact  toroid 
region  is  due  to  the  external  accelerating  field  which  is  intentionally  chosen  comparable  to 
the  internal  field.  Optimization  search  for  non-destructive  threshold  ratio  of  the  internal 
to  external  field  is  an  important  problem  for  compact  torus  accelerator  design  which  can 
be  studied  using  our  numerical  techniques. 
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List  of  Figures 


Fig.  1.  The  initial  setup  of  a  prototypical  plasma  beam  compressing  device.  Bottom 
portion:  initial  stages  of  pla.sma  beam  injection  and  compression.  Surface  compression 
shocks  are  shown. 

Fig.  2.  Compression  process  by  JxB  forces  applied  to  the  surface  of  injected  plasma 
beam.  Faster  expansion  of  plasma  at  the  front  surface  of  the  beam. 

Fig.  3.  Plasma  beam  focusing  at  the  point  of  maximum  compression  (the  ’’focal  point”). 
Top:  a  moment  before  the  focusing;  bottom:  the  focusing  instant. 

Fig.  4.  Top  portion:  the  fastest  core  part  of  the  focused  plasma  beam  leaves  the  chamber 
and  hits  a  target  located  beyond  the  computational  domain.  Bottom  portion:  formation 
of  two  consequent  temperature  maxima  near  the  beam  axis. 

Fig.  5.  A  metastable  fine  structure  of  the  beam  characterized  by  two  focal  points  and 
two  maxima  in  density  and  temperature  at  the  axis  is  formed  with  the  second  maximum 
located  near  the  outlet. 

Fig.  6.  Expansion  of  the  beam  portion  near  the  outlet  due  to  pressure  increase  induced 
by  the  beam  slowdown.  The  kinetic  energy  conversion  into  heat  causes  pressure  increase 
and  beam  expansion  which  is  most  apparent  near  the  outlet.  This  expansion  narrows 
the  channel  of  the  magnetic  flux  and  increases  magnetic  pressure  near  the  outlet  which 
accelerates  the  beam  slowdown. 

Fig.  7.  The  gap  between  the  plasma  and  the  wall  collapsed  and  magnetic  field  is  expelled 
from  the  internal  volume  of  the  chamber.  Bottom:  a  blowup  of  the  near  outlet  zone  at  the 
same  time  instant  as  shown  at  the  top.  A  strong  shock  propagates  upstream  converting 
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nearly  all  the  plasma  beam  kinetic  energy  in  the  chamber  into  heat. 

Fig.  8.  Increased  magnitude  of  the  magnetic  field  as  leading  to  premature  stagnation 
of  the  plasma  cloud  near  the  outlet.  Stagnation  is  caused  by  fast  expansion  of  the  plasma 
beam  tip  due  to  increased  rate  of  slowdown. 

Fig.  9.  A  configuration  used  to  detect  the  threshold  magnitude  of  the  magnetic  field 
which  forces  the  nearly-optimum  plasma  beam  to  focus.  Moving  plasma  initially  fills  the 
entire  chamber,  leaving  only  a  narrow  gap  for  the  magnetic  field.  The  threshold  magnetic 
field  prevents  plasma-wall  contact  which  would  otherwise  create  a  very  strong  shock. 

Fig.  10.  The  same  run:  a  moment  of  nearly  optimum  focusing.  A  bright  spot  on  the 
plasma  beam  surface  closer  to  the  inlet  is  a  trace  of  a  marginally  short  contact  between  the 
plasma  and  the  wall.  The  perturbation  caused  by  such  contact  is  still  capable  of  relaxing 
numerically;  a  very  small  relative  reduction  of  the  magnetic  field  magnitude  causes  rapid 
development  of  the  strong  shock  and  destruction  of  the  beam. 

Fig.  11.  The  same  run:  the  beam  past  the  focal  point  diverges  after  focusing. 

Fig.  12.  The  same  run:  formation  of  a  stagnated  area  near  the  outlet  shortly  after  the 
diverging  part  of  the  beam  touches  the  wall. 

Fig.  13.  The  fully  developed  structure  of  stagnated  dense  hot  low-velocity  plasma  cloud 
near  outlet.  This  moment  determines  a  ’’natural  lifetime”  of  the  magnetically  focused 
plasma  beam.  However,  the  focusing  device  efficiency  degrades  rapidly  starting  from  the 
first  implosion  of  the  shock  on  the  axis  shown  in  the  previous  Figure. 

Fig.  14.  The  plasma  beam  structure  corresponding  to  a  below-than-threshold  magnetic 
field  magnitude.  A  sequence  of  weak  shocks  collapsing  onto  the  axis  results  in  merge  and 
amplification  and,  shortly  after,  in  formation  of  shock  of  extreme  amplitude  which  forces 
the  total  stagnation. 

Fig.  15.  The  early  stage  of  the  run  with  magnetic  field  much  stronger  than  the  threshold. 

Fig.  16.  The  same  run:  a  moment  before  the  compressed  beam  leayes  the  chamber 
through  outlet.  The  first  focal  point  is  unusually  close  to  the  inlet. 

Fig.  17.  The  same  run:  a  moment  of  formation  of  the  second  focal  point  between  the 
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first  one  and  the  outlet. 

Fig.  18.  The  same  run:  formation  of  the  third  focal  point  next  to  the  outlet.  By  this 
time,  kinetic  energy  losses  in  the  beam  are  still  negligibly  small.  In  this  run  the  lifetime  of 
the  compressed  dense  beam  is  even  longer  than  in  the  runs  shown  in  Figs.  1-7  and  9-13. 

Fig.  19.  A  compact  toroid  in  prototype  accelerator  device  simulated  using  205,000 
particles,  each  one  representing  an  individual  Lagrangian  portion  of  plasma  -  a  complex 
object  which  can  develop,  change  shape  and  internal  density  and  breed  under  certain 
conditions  -  form  the  initial  body  of  the  toroid. 

Fig.  20.  A  magnified  area  around  the  toroid  shortly  after  the  accelerating  toroidal 
magnetic  flux  has  been  injected  radially  inward  through  the  inlet  above. 

Fig.  21.  The  flow/field  configuration  during  the  acceleration  stage.  Magnetic  field 
lines  leaving  the  interim  of  the  toroid  and  embedded  in  it  are  plotted  along  with  the  flow 
structures  developing  at  the  surface.  Although  the  embedded  field  is  visibly  distorted,  the 
overall  structure  is  stable. 
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